! This file is part of PanGu. All rights reserved.
module pangu_D_Alembert

    use seakeeping_kinds, only: rk
    use seakeeping_constants, only: g
    real(rk), parameter :: r13 = 1.0_rk/3.0_rk

contains

    !> 基于达朗贝尔原理，求解剖面载荷
    pure subroutine D_Alembert(mass, pressure, acceleration, xg, x, dx, width, load, acc, nstation)
        real(rk), intent(in) :: mass(nstation)  !! 质量
        real(rk), intent(in) :: pressure(nstation)  !! 外力
        real(rk), intent(in) :: acceleration(2)  !! 加速度
        real(rk), intent(in) :: xg  !! 重心坐标
        real(rk), intent(in) :: x(nstation)  !! 质量中心坐标
        real(rk), intent(in) :: dx  !! 站距
        real(rk), intent(in) :: width  !! 船宽
        real(rk), intent(out) :: load(2, 0:nstation)  !! 剖面载荷
        real(rk), intent(out) :: acc(nstation)  !! 剖面加速度
        integer, intent(in) :: nstation  !! 节点数
        integer :: i

        acc(:) = acceleration(1) - acceleration(2)*(x(:) - xg)
        load(:, 0) = 0.0_rk
        associate (force_of_inertia => mass*(g + acc))
            do i = 1, nstation
                load(1, i) = load(1, i - 1) + force_of_inertia(i) - pressure(i)*dx*width

                ! 角加速度导致的加速度(a=\alpha*r)积分存在自身内部贡献，同转动惯量计算
                load(2, i) = load(2, i - 1) + 0.5_rk*dx*(load(1, i - 1) + load(1, i)) + &
                             mass(i)*r13*(0.5_rk*dx)**2*acceleration(2)
            end do
        end associate

    end subroutine D_Alembert

    !> 求解加速度
    !> @note 采用随船坐标系
    pure subroutine get_acceleration(mass, pressure, xg, x, dx, moment_of_inertia_of_pitch, width, nstation, acceleration)
        real(rk), intent(in) :: mass(nstation)  !! 质量
        real(rk), intent(in) :: pressure(nstation)  !! 外力
        real(rk), intent(in) :: xg  !! 重心坐标
        real(rk), intent(in) :: x(nstation)  !! 质量中心坐标
        real(rk), intent(in) :: dx  !! 站距
        real(rk), intent(in) :: moment_of_inertia_of_pitch  !! 纵摇转动惯量
        real(rk), intent(in) :: width  !! 船宽
        integer, intent(in) :: nstation  !! 节点数
        real(rk), intent(out) :: acceleration(2)  !! 加速度
        real(rk) :: moment
        integer :: i

        moment = 0.0_rk
        associate (gravity => mass*g)
            acceleration(1) = (sum(pressure)*dx*width - sum(gravity))/sum(mass)
            do i = 1, nstation
                moment = moment + (gravity(i) - pressure(i)*width*dx)*x(i)
            end do
            acceleration(2) = moment/moment_of_inertia_of_pitch
        end associate

    end subroutine get_acceleration

    !> 求解中心位置
    pure subroutine get_center(mass, x, xg)
        real(rk), intent(in) :: mass(:)  !! 质量
        real(rk), intent(in) :: x(:)  !! 质量中心位置
        real(rk), intent(out) :: xg  !! 中心位置

        xg = sum(mass*x)/sum(mass)

    end subroutine get_center

    !> 求解转动惯量
    pure subroutine get_moment_of_inertia(mass, x, xg, dx, moment_of_inertia)
        real(rk), intent(in) :: mass(:)  !! 质量
        real(rk), intent(in) :: x(:)  !! 质量中心位置
        real(rk), intent(in) :: xg  !! 中心位置
        real(rk), intent(in) :: dx  !! 站距
        real(rk), intent(out) :: moment_of_inertia  !! 转动惯量

        moment_of_inertia = sum(mass*((x - xg)**2 + r13*(0.5_rk*dx)**2))

    end subroutine get_moment_of_inertia

end module pangu_D_Alembert
